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We present quantum Monte Carlo calculations of light nuclei, neutron-a scattering, and neutron 
matter using local two- and three-nucleon (3N) interactions derived from chiral effective field theory 
up to next-to-next-to-leading order (N^LO). The two undetermined 3N low-energy couplings are fit 
to the '’He binding energy and, for the first time, to the spin-orbit splitting in the neutron-a P-wave 
phase shifts. Furthermore, we investigate different choices of local 3N-operator structures and find 
that chiral interactions at N^LO are able to simultaneously reproduce the properties of A = 3,4, 5 
systems and of neutron matter, in contrast to commonly used phenomenological 3N interactions. 


Three-nucleon (3N) interactions are essential for a reli¬ 
able prediction of the properties of light nuclei and nucle¬ 
onic matter [1-5]. In quantum Monte Carlo (QMC) cal¬ 
culations phenomenological 3N interactions such as the 
Urbana [6] and Illinois [7] models have been used with 
great success [3, 8]. However, such models suffer from 
certain disadvantages: They are not based on a system¬ 
atic expansion and it was found that the Illinois forces 
tend to overbind neutron matter [9, 10]. It is therefore 
unlikely that these phenomenological models can be used 
to correctly predict the properties of heavy neutron-rich 
nuclei. 

An approach which addresses these shortcomings is 
chiral effective field theory (EFT) [2, 11-14]. Chiral EFT 
is a low-energy effective theory consistent with the sym¬ 
metries of quantum chromodynamics and provides a sys¬ 
tematic expansion for nuclear forces. It includes contri¬ 
butions from long-range pion-exchange interactions ex¬ 
plicitly and expands the short-distance interactions into 
a systematic set of contact operators accompanied by 
low-energy couplings fit to experimental data. Chiral 
EFT enables the determination of theoretical uncertain¬ 
ties and systematic order-by-order improvement; for re¬ 
cent work see Refs. [15-18]. 

Chiral EFT also predicts consistent many-body inter¬ 
actions. In Weinberg power counting, 3N forces first en¬ 
ter at next-to-next-to-leading order (N^LO) [19, 20] and 
contain three contributions: A two-pion-exchange inter¬ 
action Vc, a one-pion-exchange-contact interaction Vo, 
and a 3N contact interaction Ve- While the first is ac¬ 
companied by the couplings Ci from the pion-nucleon sec¬ 
tor, the latter two are accompanied by the couplings Cd 
and CE^ which have to be determined in A > 2 systems. 

In addition to systematic nuclear forces, reliable many- 
body methods are required to describe properties of 
light nuclei and of dense neutron matter. QMC ap¬ 
proaches, which solve the many-body Schrodinger equa¬ 
tion stochastically, are such a class of methods. Both the 
Green’s function Monte Carlo (GFMC) method and the 


auxiliary-field diffusion Monte Carlo (AFDMC) method 
rely on projection in imaginary time r, 

lim |4 't) ^ |4'o), (1) 

T—^OO 

with H the Hamiltonian of the system and |4'7’) a trial 
wave function not orthogonal to the many-body ground 
state jT'o). For a recent review of developments and 
applications of QMC methods in nuclear physics, see 
Ref. [3]. Recently, we have developed local chiral EFT 
interactions for use with QMC methods [21-24], thereby 
producing nonperturbative results for testing the chiral 
expansion scheme [22] and benchmarks for neutron mat¬ 
ter up to high density [21, 23]. However, these studies 
were limited to two-nucleon (NN) interactions only or to 
an exploratory study of neutron matter with only the 
long-range parts of the 3N interaction. 

In this Letter, we include consistent 3N interactions at 
N^LO in coordinate space [24] in GFMC calculations of 
light nuclei and n-a scattering, and in AFDMC calcu¬ 
lations of neutron matter. We fit the two couplings ce 
and Ce the ^He binding energy and low-energy n-a 
scattering P-wave phase shifts. The latter system has 
been studied using various approaches; see, for example. 
Refs. [25-27]. These observables are expected to be less 
correlated than fits to structure properties of A = 3,4 
systems because the spin-orbit and T = | components 
of the 3N interaction enter directly. 

In phenomenological 3N models, any short-range parts 
which arise from the Fourier transformation of pion ex¬ 
changes are typically absorbed into other short-distance 
structures: We retain these explicitly. We choose the 3N 
cutoff PsN = Po, where Pq is the NN cutoff, and vary the 
cutoff in the range Pq = 1.0 — 1.2 fm [21-24]. Note that 
with a finite cutoff certain ambiguities appear, including 
the specific operator form associated with the shorter- 
range interactions. In the Fourier transformation of Vd, 
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two possible operator structures arise: 


erator V on to triples with S = ^ and T = ^: 
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where Xik{r) = [5jfe(r)T(r) + cTi ■ crk]Y{r) is the 
coordinate-space pion propagator, S'ifc(r) = Scrj-rcTfe-r— 
cTi <Jk is the tensor operator, and the tensor and Yukawa 
functions are defined as T[r) = 1-1- 3/{m.„r) -|- 3/(m7r?’)^ 
and Y{r) = The smeared-out delta function 

Sr 3 n(.i^) = long-range regu¬ 
lator multiplying Y, /iong(»') = 1 — are consis¬ 

tent with the choices made in the NN interaction [21-24], 
The sum i < j < k runs over all particles 1 to A, and the 
cyclic sum runs over the cyclic permutations of a given 
triple. 

The two possible Vp structures agree in the limit of 
Rsn —>■ 0, because the delta functions then enforce i = j 
{k = j) in the first (second) term, in which case Eqs. (2a) 
and (2b) would coincide. The Vp interaction does not 
distinguish which of the two nucleons in the contact par¬ 
ticipates in the pion exchange. The second choice, Vd 2 , 
can be obtained with the exchange of a fictitious heavy 
scalar particle between the two nucleons in the contact. 
This ambiguity was also pointed out in [28]. The dif¬ 
ferences between Eqs. (2a) and (2b) are regulator effects 
and therefore higher order in the chiral expansion, but 
it is important to investigate how they affect different 
observables at this order. 

Similar effects arise in the 3N contact interaction Vp- 
Here, the main ambiguity is the choice of the 3N contact 
operator. The same Fierz-rearrangement freedom that 
allows for a selection of (mostly) local contact operators 
in the NN sector up to N^LO exists in the 3N sector at 
this order. Symmetry considerations allow the choice of 
one of the following six operators [20]: 

{l,(Ti-(Tj,Ti-Tj,(Tp(TjTpTj, 

CTi-a-jTi-Tkj [{(Ti X (Tj)-Crk][iTi X Tj)-Tk]}. 

The usual choice is ti-Tj. Here, we investigate two other 
choices: first the operator 1, and second a projector op- 


i<j k<l 

where the sums are over pairs in a given triple. In the 
infinite-momentum cutoff limit, only these S' = |,T = | 
triples would contribute to Vp due to the Pauli princi¬ 
ple. Thus, in the following we will explore three possible 
structures: 

= E (5a) 

^ i<j<k eye 

Vpi = -^^ ^ E^«3«(l'fci)<5i^3Jv(I■^J)> (5b) 

^ i<j<k eye 

Vev = H E’^'^^3«(i’fel)'^J?3iv(r*i)- (5c) 

^ i<j<k eye 

We stress that there exist other possible operator- 
structure possibilities for Vp and Vp, which will be in¬ 
vestigated in future work. 

Having specified all 3N structures, we vary the values 
of the couplings cp and cp to fit the "^He binding energy 
as shown in Fig. 1(a). We display curves for Vpi and yD 2 
using Vpr and both cutoffs i?o = 1-0 fm and i?o = 1-2 fm. 
In addition, we show curves for Vp 2 using the other two 
possible Vp structures and the cutoff i?o = 1.0 fm. For all 
of these possibilities, the stars give the values for the cou¬ 
plings which also fit P-wave n-a scattering phase shifts, 
as shown in Fig. 1(b). The resulting couplings cp and cp 
are given in Table I. In all cases (Vp) is repulsive in '^He, 
except for the case with {D2, Et) with the softer cutoff 
(i?o = 1-2 fm), where it is mildly attractive. 

For i?o = 1-0 fm and Vpr, c_d « 0 and both forms of 
Vp simultaneously fit the ^He binding energy and the P- 
wave n-a scattering phase shifts [see Fig. 1(b)]. However, 
in the softer-cutoff case Rq = 1.2 fm, Vpi and Vp 2 lead 
to different couplings. For Vpi, the splitting between 
the two P waves appears to saturate in for values of 
Cp > 2; e.g., the | phase shift for cp = 2.0, 3.0, and 5.0 
at Ecm = 1-3 MeV are each ^ 75 deg, which is ~ 35 deg 
below the i?-matrix value. Since we cannot fit the P- 
wave n-a scattering phase shifts in this case {Vpi and 
Rq = 1.2 fm), we do not consider it in the following. 
Instead, for Vp 2 and Rq = 1.2 fm, the splitting can be 
fit, as is evident from Fig. 1(b). For Vp 2 using Vpi or 
Vp-p and Rq = 1.0 fm, both the "^He binding energy and 
the P-wave n-a scattering phase shifts can be simultane¬ 
ously fit: We show only the case with Vpp in Fig. 1(b). 
There, we also show the next-to-leading order (NLO) re¬ 
sults which are a clear indication that 3N forces are neces¬ 
sary to properly describe n-a scattering. Similar results 
have been found in Refs. [29-31]. Because A = 3,4 sys¬ 
tems (futher discussed below) are largely insensitive to 
odd-parity partial waves, we find no significant depen¬ 
dence on the choice of structures in Vp. However, our 
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FIG. 1. (a) Couplings ce vs cd obtained by fitting the ^He binding energy for different 3N-operator forms. Triangles are 

obtained by using Vdi and Vet, while the other symbols are obtained for Vd 2 and three different VB-operator structures. The 
blue and green lines (lower and upper) correspond to Rq — 1.0 fm, while the red lines (central) correspond to Ro = 1.2 fm. The 
GFMC statistical errors are smaller than the symbols. The stars correspond to the values of cd and ce which simultaneously 
fit the n-a P-wave phase shifts (see Table I and the right panel). No fit to both observables can be obtained for the case with 
Po = 1.2 fm and Vdi- (b) P-wave n-a elastic scattering phase shifts compared with an P-matrix analysis of experimental data. 
Colors and symbols correspond to the left panel. We also include phase shifts calculated at NLO which clearly indicate the 
necessity of 3N interactions to fit the P-wave splitting. 


TABLE I. Fit values for the couplings cd and ce for different 
choices of 3N forces and cutoffs. 


VsN 

Ro (fm) 

Ce 

CD 

N^LO {D1,Et) 

1.0 

1.2 

-0.63 

0.0 

N^LO {D2,Et) 

1.0 

-0.63 

0.0 

1.2 

0.09 

3.5 

N^LO {D2,E1) 

1.0 

0.62 

0.5 

N^LO [D2,EV) 

1.0 

0.59 

0.0 


results in n-a P-wave scattering show a substantial sen¬ 
sitivity: Vdi appears to have a smaller effect than Vd 2 - 
In Fig. 2, we show ground-state energies and point pro¬ 
ton radii for A = 3,4 nuclei at NLO and N^LO using Vd 2 
and Vet for Rq = 1.0 fm and Rq = 1.2 fm, in compar¬ 
ison with experiment. The ground-state energies of the 
A = 3 systems compare well with experimental values. 
The ground-state energy of ^He is used in fitting cd and 
Ce, and so it is forced to match the experimental value to 
within « 0.03 MeV. The point proton radii also compare 
well with values extracted from experiment. The theo¬ 
retical uncertainty at each order is estimated through the 
expected size of higher-order contributions; see Ref. [32] 
for details. We include results from LO, NLO, and N^LO 


in the analysis using the Fermi momentum and the pion 
mass as the small scales for neutron matter (discussed 
below) and nuclei, respectively. The error bars presented 
here are comparable to those shown in Ref. [33] , although 
it is worth emphasizing that our calculations represent a 
complete estimate of the uncertainty at N^LO since we 
include 3N interactions. Other choices for 3N structures 
give similar results. 

It is noteworthy that NN and 3N interactions derived 
from chiral EFT up to N^LO have sufficient freedom such 
that n-a scattering phase shifts in Fig. 1(b) and proper¬ 
ties of light nuclei in Fig. 2 can be simultaneously de¬ 
scribed. The failures of the Urbana IX model in under¬ 
binding nuclei and underpredicting the spin-orbit split¬ 
ting in neutron-rich systems, including the n-a, system 
were among the factors motivating the addition of the 
three-pion exchange diagrams in the Illinois 3N mod¬ 
els [7]. Our results show that chiral 3N forces at N^LO, 
including the shorter-range parts in the pion exchanges, 
allow the simultaneous fit. These interactions should be 
tested further in light p-shell nuclei. 

Finally, we study the full chiral N^LO forces, includ¬ 
ing all 3N contributions, in neutron matter to extend the 
results from Ref. [24]. More specifically, we examine the 
effects of different Vd and Ve structures on the equation 
of state of neutron matter. Although these terms vanish 
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in the limit of infinite cutoff, they contribute for finite 
cutoffs. In Fig. 3 we show results for the neutron mat¬ 
ter energy per particle as a function of the density calcu¬ 
lated with the AFDMC method described in Refs. [3, 34]. 
We show the energies for Rq = 1.0 fm for the NN and 
full 3N interactions. We use Vd 2 and the three different 
Ve structures: Vet (blue band), Vei (red band), and 
Vev (green band). The error bands are determined as in 
the light nuclei case. The Vev interaction fits A = 4,5 
with a vanishing ce] hence, this choice of Ve leads to 
an equation of state identical to the equation of state 
with NN + Vc as in Ref. [24] (the projector V is zero for 
pure neutron systems), and qualitatively similar to pre¬ 
vious results using chiral interactions at N^LO [35] and 
next-to-next-to-next-to-leading order [36] . 

As discussed, the contributions of Ve and Ve are only 
regulator effects for neutrons. However, they are sizable 
and result in a larger error band. At saturation den¬ 
sity no ~ 0.16 fm“^, the difference of the central value 
of the energy per neutron after inclusion of the 3N con¬ 
tacts Vei or Vet is ^ 2 MeV, leading to a total error 
band with a range of ~ 6.5 MeV when considering differ¬ 
ent Ve structures. This relatively large uncertainty can 
be qualitatively explained when considering the following 
effects. Because the expectation value has 

a sign opposite to that of the expectation value (1) in 
"^He, ce will also have opposite signs in the two cases to 



FIG. 2. Ground-state energies and point proton radii for A = 
3,4 nuclei calculated at NLO and N^LO (with Vd 2 and Vet) 
compared with experiment. Blue (red) symbols correspond 
to Rq — 1.0 fm {Rq = 1.2 fm). The errors are obtained as 
described in the text and also include the GFMG statistical 
uncertainties. 



FIG. 3. The energy per particle in neutron matter as a 
function of density for the NN and full 3N interactions at 
N^LO with i?o = 1-0 fm. We use Vd 2 and different 3N contact 
structures: The blue band corresponds to Vet, the red band 
to Vei, and the green band to Vev- The green band coincides 
with the NN -|- 27r-exchange-only result because both Vd and 
Ve vanish in this case. The bands are calculated as described 
in the text. 


fit the binding energy. However, in neutron matter both 
operators are the same, spreading the uncertainty band. 
A similar argument was made in Ref. [37]. 

With the regulators used here, the Fierz- 
rearrangement invariance valid at infinite cutoff is 
only approximate at finite cutoff, and hence the different 
choices of Ve and Ve can lead to different results. 
The different local structures can lead to finite relative 
P-wave contributions. These can be eliminated by 
choosing Vev^ which has a projection onto even-parity 
waves (predominantly S waves). The usual nonlocal 
regulator in momentum space does not couple S and P 
waves. 

In conclusion, we find for the first time that chiral in¬ 
teractions can simultaneously fit light nuclei and low- 
energy P-wave n-a scattering and provide reasonable es¬ 
timates for the neutron matter equation of state. Other 
commonly used phenomenological 3N models do not pro¬ 
vide this capability. These chiral forces should be tested 
in light p-shell nuclei, medium-mass nuclei, and isospin- 
symmetric nuclear matter to gauge their ability to de¬ 
scribe global properties of nuclear systems. 

We also find that the ambiguities associated with 
contact-operator choices can be significant when mov¬ 
ing from light nuclei to neutron matter and possibly to 
medium-mass nuclei, where the T = | triples play a 
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more significant role. The reason for the sizable impact 
may be the regulators used here, which break the Fierz- 
rearrangement invariance, making further investigations 
of regulator choices a priority. The impact of these am¬ 
biguities in the contact operators can contribute to the 
uncertainties and needs to be studied further. 
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